arXiv:1504.01331vl [math.NA] 6 Apr 2015 


Robust split-step Fourier methods for simulating 
the propagation of ultra-short pulses in single- 
and two-mode optical communication fibers 


Ralf Deiterding and Stephen W. Poole 


Abstract Extensions of the split-step Fourier method (SSFM) for Schrodinger-type 
pulse propagation equations for simulating femto-second pulses in single- and two¬ 
mode optical communication fibers are developed and tested for Gaussian pulses. 
The core idea of the proposed numerical methods is to adopt an operator split¬ 
ting approach, in which the nonlinear sub-operator, consisting of Kerr nonlinearity, 
the self-steepening and stimulated Raman scattering terms, is reformulated using 
Madelung transformation into a quasilinear first-order system of signal intensity 
and phase. A second-order accurate upwind numerical method is derived rigorously 
for the resulting system in the single-mode case; a straightforward extension of this 
method is used to approximate the four-dimensional system resulting from the non- 
linearities of the chosen two-mode model. Benchmark SSFM computations of pro¬ 
totypical ultra-fast communication pulses in idealized single- and two-mode fibers 
with homogeneous and alternating dispersion parameters and also high nonlinear¬ 
ity demonstrate the reliable convergence behavior and robustness of the proposed 
approach. 


1 Introduction 


As computational capabilities are continuously rising, so is the demand for en¬ 
hanced networking speed. One possible approach for increasing data throughput 
is the design of networks with transmission speeds well in the Tb/s range. While the 
maximal single channel communication speed in demonstrated wavelength division 
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multiplexing systems is generally below 100Gb/s, cf. [7], we are in here concerned 
with the modeling of single- and two-mode mode optical fibers that are suitable in 
particular for long-distance data transmission. 

At present, computational models for investigating the propagation of light pulses 
in fibers have been developed primarily for pulses with a temporal half width well 
in the pico-second regime. Pulses with half widths To ^ 1 ps are sufficient for rep¬ 
resenting even on-off-key modulated bit streams with up to 100 Gb/s frequency. 
However, bit streams in the Tb/s regime can only be represented with ultra-fast 
pulses satisfying To < lOOfs. Yet, in the ultra-fast pulse regime nonlinear pulse self- 
steepening and nonlinear stimulated Raman scattering are not negligible anymore 
and an extended version of the Schrodinger-type pulse propagation equation has to 
be considered. 

Numerical solutions of the Schrodinger-type pulse propagation equation are pri¬ 
marily obtained with split-step Fourier schemes that perform spatial propagation 
steps considering firstly only the linearities in the equation by discrete Fourier trans¬ 
formation and then secondly only the nonlinear terms. While the construction of 
such split-step Fourier methods (SSFM) is very well established, cf. [1, 10], the 
topic of how to incorporate both self-steepening and Raman scattering reliably into 
the SSFM has received little attention. Here, we will describe a new class of ex¬ 
tended SSFM that properly consider the hyperbolic nature of the nonlinear sub¬ 
operator for single- and coupled two-mode optical communication fibers. 

The paper is organized as follows: In Sect. 2, we recall the governing equations 
of pulse propagation in single-mode fibers. Section 3 first discusses the construction 
principles of split-step Fourier methods and then proceeds by describing our new 
type of single-mode SSFM for ultra-fast pulses as first- and second-order accurate 
numerical schemes, cf. [5]. An ultra-fast Gaussian pulse benchmark confirming ro¬ 
bust second-order accuracy of the overall SSFM and demonstrating its application 
for simulating pulse propagation through an idealized dispersion-managed single¬ 
mode communication line are given. In Sect. 4, we describe an extended two-mode 
model for considering the simultaneous and fully coupled propagation of two ultra¬ 
fast pulses in a single fiber cable. The subsequent Sect. 5 presents a fractional step 
approach for effectively extending the derived single-mode nonlinear sub-operator 
to the corresponding system in the two-mode case. A two-mode benchmark of two 
interacting ultra-fast Gaussian communication pulses confirms the reliability of the 
method and its straightforward applicability in the dispersion-managed case is also 
shown. The conclusions are given in Sect. 6. 


2 Governing equation for ultra-fast pulses in a single-mode fiber 

The most general equation representing single-mode pulse propagation in a one¬ 
dimensional optical fiber reads 
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dA a , 




A / R{t')\A{t-t')\^dt' 


( 1 ) 

Equation (1) is derived from the electric field of the Maxwell equations, cf. [1], and 
describes the evolution of the slowly varying field envelope A{z,t) of the complex¬ 
valued signal over the propagation distance z and time t. The coefficients jSk model 
signal dispersion. Since the refractive index n of the fiber material is dependent on 
the light’s circular frequency (O, different spectral components associated to a pulse 
travel at slightly different velocities, given by c/n((o), with c denoting the speed 
of light in vacuum. This effect is mathematically modeled by expressing the mode 
propagation constant j3 in a Taylor series about the central frequency coq = 
as 

l5((o) = n{co)- = Y,^Mco-coof. ( 2 ) 

C K. 


Here, the wavelength of the injected laser light is denoted by Ao and the parameters 
a and 7 model linear signal loss and fiber nonlinearity, respectively. The function 
R{t) represents intrapulse Raman scattering, a nonlinear effect transferring energy 
from higher to lower light frequencies. Using R{t) = (1 — fR)5{t) + fRhf;{t) with 
/r = 0.18 [4] as Raman response function, applying a Taylor series expansion and 
neglecting higher order terms, Eq. (1) eventually becomes 


dz 


a ^ „ dA ^2 


li^d^A 



i d 

COo dt 


{A\A\^)-TrA 



(3) 

In general, Eq. (4) is widely accepted as a valid model for modeling the propagation 
of pulses with a half width Tq > lOfs [1]. Eor Ao = 1550nm, a typical value for the 
Raman response parameter is 7 r = 3fs. The first nonlinear term on the right-hand 
side of Eq. (3) is called the Kerr nonlinearity and the second represents nonlinear 
pulse self-steepening. 

Introducing the signal group velocity v with j3i = 1/v and using the transfor¬ 
mation T = t — zjv into retarded time T, Eq. (3) is transformed into the frame of 
reference of the pulse to read 


dA a ^ fh. d^A 


j33 d^A 

6 dT^ 


= iy{A\A\^ + iS-^{A\A\^)- TrA ^ 


(4) 


where we have also introduced 5 = (Bq *. Eor To ^ Ips, the last two terms can be 
neglected and Eq. (4) reduces to 


dA a .Pi d^A Pi d^A 
dz~^ 2 2 dT^ 6 dT^ 


iyA\A\\ 


(5) 


where j33 = 0 can be employed if Ao is not close to the zero-dispersion wavelength. 
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3 Numerical methods for ultra-fast pulses in single-mode fibers 


3.1 Split-step Fourier approach 


In order to develop a numerical solution method, Eq. (4) is commonly written in the 
form 


dz 


a 

~ 2 ~' 


.Ih , j33 


2 dT^ 6 dT^ 


A + iy(\A\^ + iS~{A\A\^)-T„^-^'^A, 






( 6 ) 


where we denote with S!{A) the operator of all terms linear in A and with •A'{A) the 
operator of all nonlinearities. Using these definitions, we write Eq. (6) in short as 


dA 

dz 


(^ + ^)A. 


(7) 


If one assumes S! and jV to be independent of z, Eq. (7) can be integrated exactly 
and the solution at z + /z reads 


A(z + /z,r) =exp(/z(^ + ,A"))A(z,7’). (8) 

The last expression forms the basis of split-step numerical methods [1]. Note, how¬ 
ever, that the operators S) and jV in general do not commute and that it corresponds 
to an 0(K) approximation to replace Eq. (8) with exp(/z^) eyzp{hjy')A{z,T). A com¬ 
monly used symmetric approximation is [24, 6] 

A{z + h,T) = exp exp(/z,A^)exp A{z,T). (9) 

Utilizing the Baker-Campbell-Hausdorff formula for expanding two non-commuting 
operators, Eq. (9) can be proven to be an 0{h^) approximation [19]. Comprehensive 
descriptions of the split-step approach for simulating pulse propagation in fibers are 
given for instance by Agrawal [1] and Hohage & Schmidt [10]. The efficiency of the 
SSEM, especially for longer propagation distances, as required for modeling optical 
communication lines, can be improved by taking solution adaptive steps in space as 
proposed by Sinkin et al. [21]. 

Alternatively, one may also construct a fractional step splitting method by solving 

r)A - 

-^ = 9A, — =^(A)A = ^(A) (10) 

dz dz 

successively, which we approximate with the symmetric fractional step method 




A{z,T), 


(11a) 



Robust split-step Fourier methods for propagation of ultra-short pulses in optical fibers 


5 


A** =A* +h.y^{A*), (lib) 

A(z + /t,7’)=exp(^^^^A”. (11c) 

Note that step (11b) is written here as a simple explicit Euler method to motivate 
the fundamental idea but schemes described below are in fact more complicated. 


3.2 Linear sub-steps 


Since the dispersion parameters j32 and (5^ are very small, discretization of the tem¬ 
poral derivatives in ^ by finite differences and approximation in physical time is 
no viable option. Instead, Fourier transformation into frequency space is commonly 
applied. The linear operator then becomes 


exp 


A{z,T)=.^ 'exp 


h (.^2 2 pi 3 



( 12 ) 


where and * denote Fourier and inverse Fourier transformation, respectively. 
In the practical implementation, discrete Fourier transformation needs to be used 
and for (O we employ the discrete frequency spectrum 


{jA(0 : j eZ A -N <j <N-1} 


(13) 


with spectral width Aco = n/{NAT). Here, it is assumed that the temporal window 
traveling with the pulse is discretized with 2N points (note that discrete Fourier 
transformation algorithms are specially efficient if the number of points is a power 
of 2), A T denotes the temporal discretization width and the temporal window has 
the extensions [—NAT, (N — l)4r]. 


3.3 Nonlinear sub-steps 

The nonlinear operator jV of the split-step method (9) is discretized in physical 
space. Utilizing |Ap = AA to eliminate 1 /A, we write aK{A) in the form 

^(A) = iy (^\A\^ + iSA^ + [iS-Tn] . (14) 

A consistent numerical method can be constructed by simply approximating the 
temporal derivatives in Eq. (14) by complex-valued first-order central differences 
and applying Eq. (9). The resulting split-step scheme would be second-order ac¬ 
curate in time and space. However, is is also clear that central finite differences 
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will result in Gibbs phenomena (cf. [13]) when strong self-steepening occurs or the 
propagation of an initially discontinuous signal needs to be simulated. 

An alternative approach for handling ^{A) is to apply forward and inverse 
Fourier transformation individually to the derivatives, cf. [16]. For instance, in (14) 

one simply replaces AdjA and with A^^^ {ico^(A)) and , 

respectively, thereby neglecting the dependence of A on T. The result is class of nu¬ 
merical operators that would generally not be consistent in the strict mathematical 
sense with ^(A) and that are not uniquely defined, with different authors arriving 
at slightly different disretizations of Eq. (14), cf. [16] and [2]. Therefore, we have 
opted to pursue a different approach, which can handle self-steepening and arbi¬ 
trary signal shapes without artificial numerical oscillations. This method is based on 
solving 


dz 


V 2 ' 2 6 dT^) 


A + /y (A|Ap + 


(A|A|2)-r«A^) 


^ yF(A) 

(15) 

within the fractional step method (11). Specific to our approach is that we discretize 
and numerically solve the complete sub-operator 


^=J^(A)= iy J^A|A|2 + iS^ (A|A|2) - TrA^^^ 


(16) 


directly. Using the Madelung transformation [17, 23] A{z,t) = \/one 
can transform Eq. (16) into the equivalent system of partial differential equations 


dz 


+ ySi 


di 

dz 


-3ySI 


dd 

^ + 7Tr 


dl 

dl 


= 0 , 

= yl 


(17a) 

(17b) 


of the real-valued quantities intensity I and phase <j). If we write the latter in the 
form 


d 

dz 



+ 

'3ySI 0 ■ 
TTr ySI 

d 

dT 

r 

= 

■ 0 ■ 


its structure as a hyperbolic advection problem 


(18) 


dz 


M(„|i 


s(q) 


(19) 


with q = (/, ^ )^ becomes apparent. The matrix M(q) has the eigenvalues Xi = 3ySI, 
Xz = ySI and a unique eigendecomposition for 1^0. Here we propose a numerical 
method for (18) that considers the characteristic information, i.e., the sign of the 
eigenvalues of M(q) for constructing one-sided (aka “upwinded”) differences for 
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the temporal derivatives, as it is required for a reliable and robust method following 
the theory of hyperbolic problems (cf. [ 22 ]). 

Again, we adopt an operator splitting technique and, instead of discretizing (19) 
directly, alternate between solving the homogeneous partial differential equation 

5,q + M(q)57-q==0 (20) 

and the ordinary differential equation 

= s(q) (21) 

successively, using the updated data from the preceding step as initial condition. 
A first-order accurate upwind scheme for (20) can be derived easily based on the 
discrete update formula [15] 

^ (M^(qi+i,q7)4q}+i/2 + M+(q,-,q,-_i)Aq"_i/2) (22) 

with Aq"^j ^2 = q/+i ~ qp where we assume a computational grid with equidistant 
mesh widths AT'm time indexed with j G Z, where —N < 7 < — 1, cf. Sect. 3.2. 

The spatial update steps are indexed by n € Nq. In general, the matrices M+ and 
indicate decompositions of M with only positive and negative eigenvalues, re¬ 
spectively. However, in the case of Eq. (18) the eigenvalues have the same sign, 
which depends solely on the sign of 7 (since / > 0). Based on (22), we construct a 
straightforward upwind scheme for Eq. (18) that reads 



(23a) 

= 0" - ^ [yTrAI] -f ySI]A(^]] , 

(23b) 


(23c) 


with 


/«^ 1 (/«+/«_i), A/j = q-iu for7> 0 , 

/« = 1 (/« , Aq = II, -1] for 7 < 0 . 

When computing the phase difference it of crucial importance to remember 
that phase is given only modulo 2n. Here, we have obtained reliable and stable 
results by ensuring that the smallest possible difference A (j)" modulo 2n is applied 
in (23b). Using the auxiliary variable 


Ae = 


0 " - <j)J_ j, for 7 > 0 , 
^j+i ~ fo*" 7 < 0- 


(24) 
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and 4 = min 10” |, |4 0” + 27r|, |4 0” — 27r| | we evaluate A 0" as 

fA9", if\Aej\=AT'l, 

A(j)’l = I A 0" + 27r, if |40” + 27r| = 4 t”, (25) 

[40;-27r, if |40”-27r| =4 t”. 

The scheme (23) is of first-order accuracy and thereby entirely free of producing 
numerical oscillations in the approximation of Eq. (15) provided that the stability 
condition 

3|7|5max{/”} A<i (26) 

is satisfied. Our present implementation guarantees (26) under all circumstances 
by having the ability to adaptively take k steps with step size Az with h — kAz 
within the central, nonlinear sub-step (lib) when required. Note, however, that for 
all computations presented in here the stability conditions (26) was always already 
satisfied for k= 1. 

To complete the algorithmic description we remark that we set /? := |A* p and 

(pj := arg(Ap after sub-step (11a) and compute A*.* = before step (11c). 

Periodic boundary conditions could be implemented by one layer of halo points. But 
note that thanks to the directional dependence, inherent to (23) and (24), it suffices 
to update only the upstream halo point, that is the one with index j = —N — 1 for 
7 > 0 and the one with j = N in case 7 < 0 before applying the upwind scheme. 


3.4 High-resolution upwind scheme 

To enable overall second-order numerical accuracy of the fractional step method 
(11), in case the solution is smooth and differentiable, it is necessary to extend the 
homogeneous nonlinear update (22) to a high-resolution scheme. For this purpose, 
we have developed a special MUSCL-type slope-limiting technique of the solution 
vector q. Originally proposed by van Leer for hyperbolic equations in conserva¬ 
tion law form [14], application to quasilinear systems is not apparent. Inspired by 
Ketcheson & LeVeque [12], we formulate our high-resolution method as 

qr' = ‘i" - ^ (M-Aq*+i/2+M+Aq*_i/2+MAq*) (27) 

with Aq*_^j/2 = q'+i - Qp “d Aq* = q^ - q'. Here, q^" 

refers to slope-limited values constructed for each component of q separately as 

qj = qj + jO’/, q^j = qj — (28) 

with reconstructed linear local slope 
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(7, = <P 






(29) 


with = qj — qj-i, ^ 7 + 1/2 = ?;+i ~ qj- latter, ^>(-) is a typical limiter 

function, where we utilize in here exclusively the van Albada limiter 

0(r) =max(0,(r2-Fr)/(l-f r^)). (30) 


To permit second-order accuracy overall, we do not utilize in (28) the discrete values 
from the previous step q" but instead intermediate values q computed as 

Qt =q"- ^ (M-4q”+i/2 + M+4q«_i/2) • (31) 

The consecutive application of (27) and (31) corresponds to an explicit 2-step 
Runge-Kutta method in the spatial update. Finally, a second-order accurate sym¬ 
metric operator splitting [24, 6] is employed to integrate Eq. (21) before and after 
the high-resolution scheme. Thanks to the simplicity of s(q) using an explicit Euler 
method for this step is equivalent to an explicit 2-step Runge-Kutta update. 

We want to point out that the first-order method (23) as well as the MUSCL-based 
second-order scheme are equally applicable for Tr—Q and especially in the singular 
case S = 0, which allows deactivation of Raman scattering and/or self-steepening if 
desired. Note that for 5 = 0 or max |/" | = 0 , the stability condition (26) is trivially 
satisfied. 


3.5 Simulation of a propagating pulse 

In order to demonstrate the described numerical method we simulate the propaga¬ 
tion of a Gaussian pulse with initial shape 

A(O,r) = y]);exp(^-^^0 (32) 

in a homogeneous fiber. The fiber is assumed to be lossless (a = 0) for simplicity 
as the omitted linear weakening of the signal is unproblematic for any numerical 
scheme. We use the SSEM in line with Eq. (11) with second-order accurate upwind- 
based nonlinear operator, cf. Sect. 3.4, and Van Albada slope-limiter (30). 

Used parameters correspond to a typical ultra-short communication pulse with 
Pq = 0.625 mW, To = 80fs, and no chirp, i.e. C = 0. The central wavelength is set to 
Xo = 1550nm, from which one computes the self-steepening parameter S = XqI 2%c, 
with c denoting the speed of light in vacuum. Raman scattering is activated with 
Tr = 3fs. Realistic fiber parameters [32 = 0.5ps^km^\ = 0.07ps^km^^ and y — 

0.1 W/m are used. Eor this configuration, the second-order dispersion length is just 
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Fig. 1 Simulated signal at 
i-max = 1 km (temporal win¬ 
dow enlarged) for Benchmark 
1. The initially Gaussian 
pulse, cf. Eq. (32), with 
half width Tq = 80 fs has 
broadened severely because 
of second-order dispersion. 
Asymmetric high-frequency 
oscillations have been added 
by third-order dispersion ef¬ 
fects. Maximal signal strength 
is reduced by a factor of 
- 13.9. 



t [ps] 


= Tq /\fi 2 \ ~ 12.8m, the third-order dispersion length is Tq /\^ 2 \ ~ 7.31m, and 
the nonlinear length is L„/ = = 16km. The pulse is assumed to travel a 

distance of just Lmax = 1 km and the simulated temporal window moving with the 
pulse has the width [—30ps,30ps —AT]. 

Figure 1 shows the computed solution using a temporal discretization of 2N 
points for = 2048 and after taking M — 100 spatial steps of equal size of h= 10m. 
Because of the very small second- and third-order dispersion lengths, typical for 
ultra-fast pulses, the final signal shape is clearly dominated by dispersion effects. 
Second-order dispersion has introduced severe pulse broadening, reducing the max¬ 
imum in power by a factor of ~ 13.9; third-order dispersion has added high- 
frequency oscillations. 

A detailed numerical analysis verifies the convergence and expected order of 
accuracy of the scheme. Starting from A = 512 and h = 4Qm (M = 25 steps), in each 
successive computation the number of Fourier modes and spatial steps is doubled. 
The numerical error at Lmax is measured for the intensity of the signal in the discrete 
maximum norm 

£oc= max \Ij-f^\jAT)\, (33) 

where a highly resolved result with N = 131,072 and M = 6400 is used as reference 
solution Figure 2 visualizes the numerical error over h and it is eminent 
that the method achieves almost perfect second-order approximation accuracy and 
reliable, robust convergence. A more detailed numerical study of the second-order 
accurate upwind-based SSFM including comparisons with several alternative nu¬ 
merical methods can be found in [5]. 


3.6 Spatially dependent fiber parameters 

Continued propagation of the pulse of Fig. 1 will invariably lead to a signal which 
has broadened to such an extent that it can not be used for digital communication. 
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Fig. 2 Numerical error over h for Bench¬ 
mark 1. The dotted line corresponds to ideal 
second order approximation accuracy. 



Fig. 3 Benchmark 2: Maximal power over 
distance with and without dispersion manage¬ 
ment. 


Yet, this problem can be compensated surprisingly easily by combining fiber sec¬ 
tions with positive and negative dispersion characteristics into a single communi¬ 
cation line. This technique is called dispersion management and has been studied 
extensively both theoretically and numerically because of its practical significance 
for long-distance fiber optical communication [18, 20, 3]. Instead of Eq. (15), one 
considers the extended variant 

dA_( a{z) Mz) , Hz) y, 

\ 2 2 dT^ 6 dT^) 

+ iy(z) [A\A\^ + iS^ (A|A|2) - (34) 

as governing equation. Adopting the practical viewpoint that the spatial numerical 
steps of any SSFM will be significantly larger than the spatial extension correspond¬ 
ing to the used temporal simulation window moving with the pulse, a straightfor¬ 
ward numerical method for Eq. (34) can be constructed by simply averaging the 
spatially dependent parameters between discrete propagation steps, i.e. by using 

Zj+Z 

i3{2,3},i = ^ J P{2,3}i^)d^ , «; = ^ j cc{^)d^ (35) 

Zj ZJ 

in the linear numerical operator (12) and by using 

Zj+h 

yj = \ J y{i)d^ (36) 


in the nonlinear operator approximating (16). 
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T [fs] Frequency 

Fig. 4 Benchmark 2: Pulse shape and spectrum after propagating 100km or experiencing 25 
soliton-like oscillations from alternating signs of dispersion parameters. 


In practice, very sophisticated dispersion management designs might be em¬ 
ployed (for instance, Guo & Huang [ 8 ] propose an exponential decrease of |j 32 | to 
accommodate better for linear loss). Here, we simply extend the example of Sect. 3.5 
and alternate the sign of ^2 and jSs every 2 km. All other parameters are unaltered 
and for an example computation we use N — 4096 and h — AQm (M = 2500) to 
simulate a pulse propagation over a distance of 100 km. In the fiber sections with 
negative dispersion parameters, the pulse deterioration is effectively reversed and 
the pulse shape mostly recovered. The pulse is undergoing a soliton-like oscillation 
with a period of 4km, which can be inferred from Fig. 3. This graphic compares the 
pulse power peak over distance in the simulation with periodic dispersion manage¬ 
ment and when the computation of the previous section is continued to a length of 
10km. In Fig. 4 are compared the shape and spectra of the initial Gaussian pulse 
and of the signal after propagating for 100 km. The observed slight signal delay and 
spectral modification is the combined effects of the nonlinearities. If 7 = 0 is used, 
the initial signal is exactly recovered. 


4 Governing equations for two interacting ultra-fast pulses 

Data throughput can be increased significantly if multiple optical fields of differ¬ 
ent wavelengths propagate simultaneously inside the fiber. However, these fields 
would interact with one another through all the fiber nonlinearities. Additionally 
if three or more fields are initially present, even new signal fields can be induced 
(aka four-wave mixing [1]). Therefore, we consider in the following only the case 
of two interacting signal fields propagating through an optical fiber, for which there 
is already some agreement about the structure of the governing equations in the lit¬ 
erature [11]. Extensions of the ultra-fast pulse propagation equation (3) to three or 
more interacting fields are still a topic of active research. 

We assume two pulses at carrier frequencies ^nd two nonlinear con¬ 

stants 7 i, 72 . It is further assumed that the cross-phase modulation of each frequency 
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can be expressed for all higher order nonlinear terms by positive factors Bi, B 2 , the 
cross-phase modulation in the Kerr nonlinearity by factors Ci, C 2 . Extending Eq. (3) 
accordingly, we use the model equations 


dA\ 

dz 


dA2 

dz 


ar d^Ai n . i2 ^ 1 . i2\ . 

—^"^^6 +C 1 IA 2 I )Ai 


CO, 


7i 

( 1 ) 


5(|A 


iPAi 


Bi 


^(|A2pAi) 


dt^ 
-iJiTr 


dt ' dt 

dA2 .PP d^A2 , PP d^A 


d\A 


'2 5|A2|21 


dt 


ABv 


dt 


At, 

(37a) 




dfi 


+ ' 


CO, 


72 

( 2 ) 


5(|A2pA2) 5(|AipA2) 

-f £>2 


dt 


dt 


6 dt^ 

1 

— iJiTr 


+ iy 2 (IA2P-1-C2IA1 p) A2 


^|A2p , „ 5|AiP^ 


dt 


AB 2 - 


dt 


A 2 . 

(37b) 


Note that (37) encompasses the model actually adopted for simulation by Kalithasan 
et al. in [11]. Using pP = l/vj and the transformation T = t — pPz into retarded 
time yields 


dAi 

dz 


dA2 

dz 


^A 

2 

71^1 


.Pp'’ d^Ai pP ^3^1 


5(|A: 


pA 


dT 


2 dT^ 


dT 


dT^ 
-iYiTr 


+ *7i (|Aip-l-Ci|A2p)Ai 


^|AiP 

dT 


Bi 


dlM 

dT 


21 


At, 

(38a) 


«2. .^A2 ^PP d^A2 PP d^A2 


2 A 2 5 t 2 


dT^ 


Y 2 S 2 


5(|A2pA2) 5(|AipA2) 

+ B 2 


dT 


dT 


6 dT^ 

— iYiTr 


-1-/72 (|A2p-|-C2|Aip)A2 

^|A2p 


dT 


+ B 2 


5|Aipl 


dT 


M, 

(38b) 


with 5 = (vi — V 2 )/(viV 2 ) representing the group velocity mismatch between both 
fields. As before we use 5^ = 1 /CoP fovk= 1,2. 

In the regime of pico-second pulses, that is for pulses with Tq ^ 1 ps, two-mode 
extensions of Eq. (5) are well established. Setting Sk = 0, Tr — 0 and using Q = 2 
in (38), we obtain the frequently used [1] cross-phase modulation model 

Aid- /71 (|Aip-|-2|A2p)Ai, 


5Ai^/_ai _.pP^ pP^\ 
2 ^ 2 dT^ 6 dT^ j 


(39a) 
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dA2 _( az d .pP 
dz 2 ^ dT ^ 2 dT^ 

'■ -V- 

^{2) 




A2 + 


m(|A2p + 2|Aip)A2, 

^(2) 

(39b) 


which we write as 

^ +^(i)(Ai,A 2)) Ai, ^ +.A"(^((Ai,A 2)) Az. (40) 


5 Numerical methods for two interacting ultra-fast pulses 
5.1 Extended split-step Fourier method 

Taking advantage of the fact that the linear operators only need to be applied 
to each field Aj^, a SSFM for approximating solutions of system (39) - in line with 
Eq. (9) - is easily constructed as 

Al =exp Ai, Aj = exp Az, (41a) 

Af =exp(^/i^(^)(A|,A^))A|, A^* = exp (^/i^P((Ar,A^)) AJ (41b) 

Ai(z + /i) =exp Aj*, A 2 (z + /i) = exp A”. (41c) 

Obviously, the numerical operators of (41b) and (41c) acting on each fields can be 
executed consecutively. The linear operator is identical to (12). For we 
have 


exp 




exp 


^ ^ -^2^^ 2 3 «2 


^Az. 


(42) 

A second-order accurate scheme can be expected if (41b) is replaced with a sym¬ 
metric splitting scheme such as 


At = expQ^(')(Al,A^)^Al, (43a) 

Ar = exp {h.A^^^\A\,A* 2 )^ A|, (43b) 

Ar=expQ^(')(At,Ar))A|. 


(43c) 
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5.2 Nonlinear sub-steps 


While the derivation of a SSFM for the simplified system (39) is apparently 
a straightforward task, formulation of a reliable numerical method for the sys¬ 
tem of propagation equations for two coupled ultra-fast pulsed signals, (38), is 
more involved. In particular, when the equations of (38) are written in the form 
one quickly finds that due to the cross-phase coupling the 
factor 1 /Aj^ of the self-steepening term cannot be eliminated from as it was 
done to obtain Eq. (14). This leaves a singularity in the operator for vanishing sig¬ 
nals and neither the centered difference method nor particularly an ad hoc Fourier 
transformation technique, sketched both in the beginning of Sect. 3.3, are available 
anymore for numerical method construction. However, we will demonstrate subse¬ 
quently how our upwind-based discretization technique of Sect. 3.3 can be easily 
extended to (38), yielding a reliable and robust numerical method. 

We start the derivation of the method by inserting the linear operators from (39) 
into (38) to obtain 


dA 

dz 


y = +Ck\Ai\^')Ak 


-JkSk 


d{]A,\^A,) d{\Ai\^A,) 

+ i>k- 


dT 


dT 


- ijkTR 


d\Ak\^ 

dT 


+ Bk 


d\Ai\ 
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dT 


Ak 


(44) 


for k,l € {1,2} and k 7 ^ /. In analogy to Sect. 3.3, we assume a fractional step 
approach in the spirit of Eq. (11) that considers the linear operators with the update 
steps (41a) and (41c) and approximates the nonlinear sub-operator equations 


dAk 

dz 


^^’^\Ak,Ai) = ijk {]Ak\^ + Ck\Ai\^)Ak 


-JkSk 


d{]Ak?Ak) d{\A,\^Ak) 

+ Bk- 


dT 


dT 


- ijkTR 


d\Ak\^ 


dT 


dT 


Ak- 


_(45) 

Using again Madelung transformation for each field, i.e. Ak{z,t) = \/tk{z,t)e“i‘A^’‘\ 
we obtain the transport equations for the intensities Ik and the phases (j)k instead of 
(45) as 


d^k 

dz 


^JkSkih+Bkli)^ 


= 0 , 

= Jk {tk + Ckh) ■ 


(46a) 

(46b) 


The latter defines a single system of advection equations that couples the fields Ak 
and A/. Using the state vector u = (/i, 0i ,/ 2 , ^ 2 )^, this system reads 
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with matrix 


B(u) = 


>151(3/1+B1/2) 0 27i5iBi/i 0 

JiTr 'Y\S\{h+B\l2) JiTrBi 0 

^72S2B2h 0 72S2{^h+B 2 I 1 ) 0 

72TRB2 0 72TR 72S2ih + B2h) 


and right hand side 

r(u) = (0, 71 (/i + C 1 / 2 ), 0, 72 {h + C 2/1 ))^. (49) 

In order to verify the hyperbolicity of Eq. (47) and for constructing an upwind 
scheme, one would require the eigendecomposition B = RAR^^ However, the nec¬ 
essary linear algebra is very involved and is additionally complicated by the singular 
cases 4 = 0 , which have to be considered separately in order to construct a gener¬ 
ally robust numerical scheme. To simplify the latter, we have opted to use a splitting 
approach and update the fields Aj^ and A/ successively. Instead of solving the com¬ 
bined system (47) we construct an approximation to (46) under the assumption that 
7/ is independent of z. Proceeding then as in Sect. 3.3, we write (46) as the advection 
system 

^ ^ (qC^)), (50) 

with vector of state q^*^^ = ( 4 , ^ji,//)^, matrix 

7kSki^h+BkIi) 0 27kSkBkIk 

M(k)(q(k))= 7kSk{h+BkIi) 7kTRBk (51) 

0 0 0 


and source term 

s(k)(q(k)) ^ (0,7>4 + C,/i),0)^. (52) 

The non-zero eigenvalues of M*^k) ^re 7kSk{'iIk+ Bkh) and 7kSk{Ik+ Bkh). Since 
Bk > 0 and 4 /; ^ 0 hold true, both eigenvalues have again the same sign, solely 
determined by the sign of 7k. Following the upwind approach again we construct a 
first-order accurate method for (50) as 

4”j ‘ = 4"; - -^7kSk [(3>,. +/?,/;>A4> + 2B,»/,”j , (53a) 

= K.j-^7k [r>A4>. + /?,A/>)+5>>^.+B>;;^.)44>] , (53b) 

Ky=Ky + h 7 k {iiy +q/>) , ( 53 c) 

where 

Jn _ ^ (jn .jn \ ath _ jti jn 

k/i,j - 2 yk/ij +h/ij-\) > ^4//,; - h/ij-h/i,j-i 


for 7k > 0 , 
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jn I jn 
h/ij+^k/ij+i 


A Tit — jn _jn 


for Yk < 0 . 


As before, is evaluated modulo In using Eqs. (24) and (25) and the stability 
condition reads 

I Yk I Stmax { 34 j + Bklij } ^ ^ 1 • (54) 

By construction, the single-field upwind method (53) computes only new values for 
4 and while the intensity of the other held, ![, is assumed to remain unchanged. 
In order to achieve an update of both helds, and thereby approximation of (47), we 
apply the single-held upwind scheme within a symmetric fractional-step splitting 
method, i.e. 


Al + ’^^m{Al,A*2), 

(55a) 

A*2+hJ^^^HA\,A*2), 

(55b) 

A* + ^^W(A\,A*2*). 

(55c) 


A symmetric SSFM is obtained by applying expressions (41a), (55) and (41c) after 
another. Finally, the high-resolution technique, described in Sect. 3.4, is adopted 
to implement a second-order accurate approximation to where we presently 

apply slope-limited reconstruction to 4 and but not to /;. 


5.3 Simulation of two interacting propagating pulses 

We use a conhguration with very strong nonlinearity and thereby nonlinear pulse 
interaction to assess the reliability of the derived two-mode method. A hber without 
linear loss and third-order dispersion is assumed, i.e. ap 2 =0 and jSj ’ ^ = 0 , and 
Raman scattering is also deactivated by setting Tr — 0. To enforce a strong inhu- 
ence of the nonlinearities we use = 4 x 10^^ps^km^\ 7 i = 1W/m, and 72 = 
1.2W/m. Two unchirped pulses in the range of ultra-short communication pulses 
with = 80fs and power levels of = 0.625mW and = 0.3125mW 

are used. The hrst central wavelength is set to = 1550nm and the second to 

( 2 ) 

Aq ' = 1300nm. The group velocity mismatch parameter is set to 5 = 0.015625 fs/m 
and the cross-phase modulation parameters read Pi 2 = C 72 = 2 . 

For this configuration, the second-order dispersion length is = 160 km and the 
nonlinear lengths = 1. 6 km and = 2.667km, respectively. The approximate 

optical shock distances [1], = •\/eL|jJ’^((BQ*’^(ro/(3-\/2) are ~ 60.491km and 

~ 120.216km, respectively. We use a propagation distance of Lmax = 64km, yield¬ 
ing a temporal shift of the second pulse by exactly 1 ps, and the temporal window 
has the width [—4ps,4ps — AP], 
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Fig. 5 Benchmark 3: Pulse shape and spectrum after propagating 64 km of two individual highly 
nonlinear ultra-short single-mode pulses (solid lines) and when the two pulses are interacting with 

one another in a two-mode fiber. The upper row corresponds to Pulse 1 with = 0.625mW; the 

( 2 ) 

lower row to Pulse 2 with Pq ' = 0.3125mW and S = 0.015625fs/m causing the pulse to arrive 
1 ps earlier. 


In Fig. 5 is shown the computed solution using a temporal discretization of 
2N points for N = 2048 and after taking M — 3200 spatial steps of equal size of 
h = 20 m. Additionally are shown the solutions if each pulse travels individually. 
These solutions are computed by keeping all other parameters unchanged while 
setting = 0 and = 0, respectively. If only a single field is present, our two¬ 
mode SSFM is identical to the previously developed single-mode SSFM, which was 
confirmed to be second-order accurate in Sect. 3.5. Note that the single-mode so¬ 
lution of Pulse 1 was also used as a detailed computational benchmark in [5] and 
is thereby available as a reference. From Fig. 5 it can be seen that the two non¬ 
interacting single-mode pulses exhibit a very similar shape and spectrum. However, 
in the two-mode model particularly the faster and weaker second pulse is signifi¬ 
cantly altered. Pulse 2, visualized in the lower row of Fig. 5, experiences consid¬ 
erable signal steepening from cross-phase modulation, which can be inferred espe¬ 
cially from its spectrum. 

We use the same technical approach as in Sect. 3.5 to quantify the numerical error 
and order of accuracy of the two-mode SSFM. We double the temporal resolution 
consecutively starting from A = 512uptoA = 16,384 and simultaneously divide 
the spatial step size by a factor of 2 respectively, starting with /z = 80 m (M = 800 
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Fig. 6 Numerical error Eoo 
over h for Benchmark 3. The 
respective error of Pulse 1 
is marked with solid lines, 
the respective error of Pulse 
2 with dotted lines. Single¬ 
pulse simulation results are 
indicated with open squares, 
the fully coupled simulations 
are marked with closed cir¬ 
cles. The upper broken line 
corresponds to an order of 
accuracy ^ 1.47, the lower 
one to an order of accuracy of 
- 2 . 20 . 



steps). The numerical error at Lmax is measured for /(j 2 ) in the maximum norm, cf. 
Eq. (33), where results computed with N = 32,768 and h= 1.25 m (M = 51,200) 
are used as respective reference solutions. The computational errors of a series of 
fully coupled two-mode results as well as the errors of single-mode computations 
(cf. Fig. 5) of both individual pulses are plotted in Fig. 6 . In general, the example 
confirms that the proposed two-mode SSFM converges reliably and robustly even 
for a highly nonlinear coupled problem and performs identical beside round-off 
errors to the single-mode method of Sect. 3.4 for uncoupled individual pulses. While 
the single-mode SSFM with limiter (30) actually achieves slight super-convergence 
in this test case (the measured order of accuracy is ~ 2.20), the two-mode SSFM of 
Sect. 5.2 with same limiter yields an approximate order of accuracy of ^ 1.47. One 
might attribute this behavior to the fractional step splitting treatment of the nonlinear 
operator, (55), however increasing the number of spatial steps up to a factor of 8 to 
possible reduce the splitting error of the nonlinear sub-operator resulted only in 
marginally smaller numerical errors for this test case. 


5.4 Spatially dependent fiber parameters 

As final test case, the coupled propagation of the two Gaussian pulses of the previous 
benchmark through the dispersion-managed communication line of Sect. 3.6 is con¬ 
sidered. Fike in Sect. 3.6 we assume an optical communication line of 100km length 
with dispersion parameters = 0.5ps^km^* and = 0.07ps^km^', 

which all change sign every 2km, and o:(i 2 ) = 0, 7(1 2 = 0.1 W/m, and Tr = 3fs. As 

before, the parameters of the two unchirped Gaussian pulses are = 0.625 mW, 

= 0.3125mW, and = 80fs. The wavelengths are again = 1550nm 
and = 1300nm. The group velocity mismatch is 5 = 0.015625fs/m and cross¬ 
phase modulation parameters are B 12 = 6 ^ 1,2 = 2. The same computational param- 
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Fig. 7 Pulse shape and spectrum of two coupled pulses after propagating 100 km through the 
idealized dispersion-managed fiber of Benchmark 2. 


eters are used as in Sect. 3.6: The temporal window has the width [—30ps,30ps — 
A T] and N = 4096, h = 40 m are applied. 

During propagation both pulses are experiencing almost undisturbed soliton- 
like oscillations every 4 km. Figure 7 compares the final signal shapes and spectra 
with the respective initial ones, where Pulse 2 has been shifted for visualization by 
—1.5625ps. Both pulses are delayed by roughly lOfs but the signal shape is quite 
well preserved; the spectral alteration being rather moderate in both cases. In the 
left graphic of Fig. 7 Pulse 1 and 2 are easily distinguished; in the right graphic the 
final spectra of Pulse 1 and 2 are specially indicated. 

Finally, we comment on typical run times of the proposed split-step Fourier 
methods. Our implementation is in FORTRAN 90 and uses the Netlib NAPACK 
Fast Fourier Transformation (FFT) routines, which are coded in FORTRAN 77 [9]. 
Compiled with usual optimizations, the two-mode computation of Fig. 5 required 
~ 48 seconds on a single Intel Xeon E5 CPU with 2.1GHz. Dependence on the 
number of Fourier modes N as well as the number of spatial steps M is linear and 
each computation of the convergence analysis of Fig. 6 is therefore four times more 
expensive than the next coarser one. On the same CPU, the dispersion-managed 
two-mode simulation of Fig. 7 ran for ~ 100 seconds, its single-mode analogue 
of Fig. 4 required ^ 40 seconds. These moderate run times and the given results 
provide evidence for the relevance of the proposed numerical methods for practical 
long-distance fiber optical communication line design. 


6 Conclusions 

Reliable extensions of the classical SSFM into the regime of ultra-fast pulses have 
been derived and demonstrated for typical Gaussian communication pulses in highly 
nonlinear and dispersion-managed long-distance optical fibers. The primary diffi¬ 
culty in this regime lies in the appropriate mathematical treatment of the additional 
nonlinear terms modeling signal self-steepening and stimulated Raman scattering. 
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For the case of the single-mode equation (3) and the two-mode system (37) it was 
shown that under Madelung transformation all nonlinearities can be effectively com¬ 
bined into an inhomogeneous system of advection equations of the signal intensi¬ 
ties and phases. Following upwind and slope-limiting ideas, originally developed 
in the context of supersonic hydrodynamics, a robust numerical method is then de¬ 
rived for the single-mode nonlinear sub-operator and incorporated into a symmetric 
SSFM. Reliable convergence and numerical approximation accuracy of second or¬ 
der is demonstrated for the overall method. While it would be principally feasible 
to apply the exact same approach to the two-mode case and the correspondingly 
derived four-dimensional system (47), we have opted for now for a mathemati¬ 
cally less involved fractional step approach and apply two single-field nonlinear 
sub-operators successively to approximate the solution of (47). This single-field 
sub-operator is derived as a straightforward extension of the slope-limited upwind 
method for the single-mode case. Incorporated into a two-mode SSFM, the overall 
numerical scheme converges reliably, yet, in a highly nonlinear test case only an 
order of accuracy of ~ 1.5 is obtained. Future work will concentrate on develop¬ 
ing an unsplit scheme for (47). It is expected that such a method should obtain an 
order of accuracy close to 2 while being of comparable computational expense and 
robustness as the two-mode SSFM proposed in here. 
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